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Thermosphere  Density  Variability,  Drag  Coefficients,  and  Precision  Satellite  Orbits 


Craig  A.  McLaughlin,  Dhaval  Mysore  Krishna,  Piyush  M.  Mehta,  Travis  Lechtenberg,  Andrew  Hiatt, 

Eric  Fattig,  and  Travis  Locke 

Abstract 

Precision  orbit  ephemerides  (POE)  are  used  to  estimate  atmospheric  density  along  the  orbits  of  CHAMP 
(Challenging  Minisatellite  Payload)  and  GRACE  (Gravity  Recovery  and  Climate  Experiment).  The 
densities  are  calibrated  against  accelerometer  derived  densities  and  considering  ballistic  coefficient 
estimation  results.  The  14-hour  density  solutions  are  stitched  together  using  a  linear  weighted  blending 
technique  to  obtain  continuous  solutions  over  the  entire  mission  life  of  CHAMP  and  through  2011  for 
GRACE.  POE  derived  densities  outperform  the  High  Accuracy  Satellite  Drag  Model  (HASDM), 
NRLMSISE-00  model,  and  Jacchia  71  model  densities  when  comparing  cross  correlation  and  RMS  with 
accelerometer  derived  densities.  Cross  correlations  for  all  other  data  sets  with  accelerometer  derived 
densities  are  lower  when  the  satellite’s  orbit  planes  are  near  the  terminator.  Additional  research  showed 
that  the  high  frequency  variations  in  density  observed  only  be  accelerometers  have  little  effect  on  orbit 
propagation  accuracy.  CHAMP  and  GRACE  precision  orbit  data  were  decimated  to  8  and  15  minutes  per 
orbit  and  had  noise  levels  up  to  100  m  added  to  them  to  observe  the  effects  of  limited  data  of  lower 
accuracy  on  density  estimation.  As  expected,  the  accuracy  decreased  as  higher  levels  of  noise  were  added 
to  the  data,  but  the  even  with  limited  data  of  moderate  accuracy,  the  estimated  densities  were  of  usable 
accuracy.  The  Direct  Simulation  Monte  Carlo  was  used  to  estimate  drag  coefficients  for  the  GRACE, 
Stella,  and  Starlette  satellites.  The  resulting  drag  coefficients  showed  good  agreement  with  existing  drag 
coefficient  techniques.  The  models  were  also  used  to  develop  a  new  empirical  drag  coefficient  model  for 
GRACE,  which  can  be  easily  extended  to  other  satellites.  The  empirical  model  allows  computationally 
efficient  determination  of  drag  coefficients  for  complex  satellites. 

1  Introduction 

Atmospheric  density  and  drag  coefficient  modeling  have  long  been  among  the  greatest  uncertainties  in  the 
ability  to  predict  the  motion  of  low  Earth  orbit  (LEO)  satellites,  which  is  a  critical  aspect  of  space 
situational  awareness  (SSA).  Accurate  density  and  drag  coefficient  calculations  are  required  to  provide 
meaningful  estimates  of  the  atmospheric  drag  perturbing  satellite  motion.  These  effects  increase  with 
lower  altitude  orbits  and  also  with  higher  effective  area  and  lower  mass  satellites.  The  proposed  effort 
will  use  precision  satellite  orbits  to  examine  thermospheric  density  changes  and  examine  techniques  for 
improving  satellite  drag  coefficient  modeling. 

The  long  term  objectives  of  the  PI  are  to  improve  capabilities  in  the  area  of  SSA  by  improving  orbit 
estimation  and  prediction  for  LEO  satellites,  better  understanding  the  density  fluctuations  in  the 
thermosphere,  and  improving  the  transfer  of  knowledge  between  the  aeronomy  and  orbital  mechanics 
communities.  The  specific  objectives  of  this  proposed  effort  were  to: 

1 .  Better  understand  the  density  variations  in  the  thermosphere  by  using  precision  orbit  ephemeris  (POE) 
derived  density  data  from  several  satellites  to  simultaneously  examine  the  fluctuations  in  atmospheric 
density  over  time  scales  ranging  from  two  hours  to  one  day  and  correlate  these  fluctuations  with  specific 
upper  atmospheric  phenomena.  How  does  atmospheric  density  change  in  both  time  and  in  space?  How 
do  these  changes  affect  the  motion  of  satellites?  These  are  some  of  the  fundamental  questions  that  this 
research  will  address. 

2.  Better  understand  the  relationship  between  drag  coefficient  models  and  orbit  determination  and 
prediction  accuracy.  This  objective  will  examine  the  limitation  of  existing  theories  for  satellite  drag 


coefficients.  In  addition,  the  research  will  examine  whether  modern  computational  fluids  techniques  are 
applicable  to  this  problem  and  worth  the  increased  computational  burden  they  impose. 


2  Density  Estimation 
2.1  Background 

McLaughlin  (2005)  gives  an  introduction  to  the  neutral  atmosphere  and  the  time  varying  effects  on  the 
thermospheric  and  exospheric  density.  These  time  varying  effects  include  solar  rotation,  solar  cycle, 
diurnal  variations,  magnetic  storms  and  substorms,  gravity  waves,  winds  and  tides,  and  long-term  climate 
change.  Vallado  (2007)  gives  an  introduction  to  the  basic  variations  in  density  as  well  as  the  most 
commonly  used  density  models  in  orbit  determination.  Sabol  and  Luu  (2002)  give  a  summary  of  the 
drivers  of  atmospheric  density  variations  and  discuss  some  of  the  difficulties  connected  with  the  temporal 
resolution  of  various  proxies  utilized  by  the  empirical  density  models.  Marcos  et  al.  (2003)  supply  an 
overview  of  research  addressing  the  inaccuracies  in  modeling  satellite  drag. 

Two  primary  categories  of  research  exist  to  address  the  problems  of  modeling  atmospheric  density  for 
satellite  drag.  The  first  category  is  dynamic  calibration  of  the  atmosphere  (DCA)  and  the  second  is  using 
accelerometers  onboard  satellites  to  measure  non-conservative  accelerations,  which  includes  drag.  DCA 
utilizes  the  observed  motions  of  a  large  number  of  satellites  in  order  to  estimate  large-scale  density 
corrections  to  an  existing  atmospheric  density  model  (Bowman  et  al.,  2007;  Bowman,  2004;  Cefola  et  al., 
2003;  Storz  et  al.,  2005;  Wilkins  et  al.,  2007a,  2007b;  Yurasov  et  al.,  2004;  Yurasov  et  al.,  2008). 
Dynamic  calibration  of  the  atmosphere  provides  a  significant  improvement  to  empirical  density  models 
but  with  several  disadvantages.  First,  a  DCA  approach  is  designed  to  run  internal  to  a  particular  orbit 
determination  scheme  with  the  resulting  atmospheric  density  corrections  only  applying  to  a  specific  time 
period.  Therefore,  those  using  a  different  orbit  determination  scheme  must  rely  on  that  particular  system 
for  updates  to  atmospheric  density  corrections  as  well  as  requiring  a  complete  archive  of  density 
corrections  for  a  given  problem.  The  second  limitation  of  DCA  approaches  is  the  limited  spatial  and 
temporal  resolution  of  the  atmospheric  density  corrections.  The  corrections  allow  the  baseline  density 
model  to  characterize  atmospheric  density  variations  in  terms  of  several  hours  to  days  but  not  in  shorter 
time  scales.  A  temporal  limitation  is  introduced  by  the  use  of  a  daily  solar  flux  and  averaged  3-hour 
geomagnetic  indices  as  input  values  into  a  DCA  scheme.  The  use  of  this  input  data  does  not  permit  the 
baseline  atmospheric  density  model  to  properly  represent  variations  that  occur  within  the  averaging 
interval  of  these  indices.  Further  difficulty  arises  in  DCA  approaches  because  of  the  predominate  use  of 
two-line  element  sets  of  a  large  number  of  low  Earth  orbit  objects  as  observational  inputs  for  a  given 
DCA  approach.  Reliance  on  two-line  elements  results  in  reduced  accuracy  density  corrections  with 
restricted  temporal  resolution.  The  High  Accuracy  Satellite  Drag  Model  (HASDM)  (Storz  et  al.,  2005) 
utilizes  radar  observations  of  low  Earth  orbit  objects,  but  the  accuracy  of  radar  observations  is  still  lower 
than  those  obtained  by  precision  orbit  ephemerides  (POE)  or  satellite  laser  ranging  (SLR).  In  addition, 
the  radar  observation  data  are  not  generally  available. 

The  second  category  of  research  for  improving  atmospheric  density  knowledge  is  using  accelerometers 
onboard  satellites  to  measure  non-conservative  accelerations,  which  can  be  utilized  to  estimate  density. 
The  accelerometer  data  allows  the  separation  of  gravitational  forces  from  non-conservative  forces 
including  Earth  radiation  pressure,  solar  radiation  pressure,  and  drag.  Use  of  accurate  radiation  pressure 
models  permits  the  drag  acceleration  and  resulting  estimated  density  to  be  accurately  calculated  with 
precise  temporal  resolution.  The  accelerometer  data  are  extremely  precise  but  only  available  for  a  few 
satellites.  This  represents  an  extreme  opposite  in  terms  of  accuracy  and  data  availability  compared  with 
the  use  of  two-line  element  sets.  The  availability  of  accelerometer  measurements  is  currently  limited  to 
the  Challenging  Mini-Satellite  Payload  (CHAMP)  and  the  two  Gravity  Recovery  and  Climate  Experiment 
(GRACE)  satellites.  Other  satellites  with  accelerometers  have  flown  in  the  past.  An  example  is  given  in 


Rhoden  et  al.  (2000),  in  which  the  accelerometer  data  from  the  Satellite  Electrostatic  Triaxial 
Accelerometer  (SETA)  experiment  flown  at  an  altitude  of  about  200  km  was  used  to  examine 
thermospheric  density  and  confirmed  that  energy  from  magnetic  storms  deposited  at  high  latitudes  created 
a  “density  bulge”  that  propagated  toward  the  equator  and  then  toward  the  opposite  poles. 

Konig  and  Neumayer  (2003)  and  Bruinsma  and  Biancale  (2003a,  2003b)  published  some  of  the  early 
results  for  estimating  density  using  CHAMP  accelerometer  data  with  additional  atmospheric  density 
values  derived  using  CHAMP  accelerometer  data  by  Bruinsma  et  al.  (2004)  and  Nerem  et  al.  (2003). 
Schlegel  et  al.  (2005)  examined  polar  region  density  structures  by  using  CHAMP  accelerometer  data. 

Joint  work  between  researchers  at  the  University  of  Colorado  and  the  Centre  National  d’Etudes  Spatiales 
(CNES)  utilized  CHAMP  and  GRACE  accelerometer  data  to  examine  density  variations  created  during 
solar  and  geomagnetic  events  (Bruinsma  et  al.,  2006;  Bruinsma  and  Forbes,  2007;  Forbes  et  al.,  2005; 
Sutton  et  al.,  2005;  Sutton  et  al.,  2006;  Sutton  et  al.,  2007).  Tapley  et  al.  (2007)  describe  a  technique  to 
estimate  atmospheric  density  from  the  GRACE  accelerometer  data.  The  CHAMP  and  GRACE  satellites 
have  contributed  vast  amounts  of  information  regarding  the  upper  atmosphere  and  are  capable  of  adding 
even  more.  Unfortunately,  these  three  satellites  only  provide  limited  spatial  coverage  at  any  given  time 
and  only  at  low  altitudes. 

The  research  presented  in  this  report  is  another  step  toward  the  goal  of  combining  accurate  data  with  good 
spatial  coverage  obtained  from  a  large  number  of  satellites.  Many  satellites  currently  possess  Global 
Positioning  System  (GPS)  receivers  that  when  combined  with  precision  orbit  determination  provide 
position  accuracies  of  a  few  centimeters  to  within  a  few  meters.  This  research  estimates  atmospheric 
density  by  using  the  precision  orbit  data  of  these  satellites  in  a  precision  orbit  determination  scheme. 
CHAMP  and  GRACE  POEs  are  used  to  estimate  atmospheric  density.  Use  of  POE  data  results  in 
increased  accuracy  from  a  smaller  number  of  satellites  as  compared  with  two-line  element  sets.  When 
compared  with  accelerometer  data,  POE  data  provide  a  larger  number  of  available  satellites  with  reduced 
accuracy.  Several  papers  have  examined  the  use  of  GPS  receiver  or  SLR  observations  for  estimating  non¬ 
conservative  accelerations.  One  such  approach  is  given  by  Doornbos  et  al.  (2005),  where  a  type  of 
differential  correction  was  examined  by  using  two-line  element  sets  in  a  traditional  DCA  scheme  along 
with  a  small  number  of  satellites  with  precision  orbit  data.  Another  approach  described  as  GPS 
accelerometry,  utilizes  GPS  receiver  data  to  estimate  non-conservative  forces  as  empirical  accelerations 
(van  den  Ijssel  et  al.,  2005;  van  den  Ijssel  and  Visser,  2005,  2007).  Through  this  method  the  in-track  and 
cross-track  accelerations  derived  from  the  CHAMP  accelerometer  may  be  reasonably  determined  with  a 
temporal  resolution  of  20  minutes  or  greater.  Montenbruck  et  al.  (2005)  utilized  batch  and  Kalman  filter 
estimation  techniques  in  order  to  investigate  the  reconstruction  of  empirical  accelerations  of  the  GRACE- 
B  satellite.  Both  the  batch  and  Kalman  techniques  demonstrated  similar  overall  variations  in  the 
empirical  accelerations,  but  with  a  scale  difference  between  the  acceleration  magnitudes  obtained  from 
the  two  methods.  Willis  et  al.  (2005)  used  Doppler  Orbitography  and  Radiopositioning  Integrated  by 
Satellite  (DORIS)  and  SLR  data  to  examine  density  variations  in  the  800-900  km  and  1300-1400  km 
ranges  in  the  thermosphere  during  periods  of  enhanced  geomagnetic  activity.  The  study  found  significant 
errors  in  the  atmospheric  models,  but  these  errors  were  greatly  reduced  with  additional  enhanced  data 
processing.  DORIS  is  yet  another  way  of  obtaining  highly  accurate  satellite  state  vectors,  and  allows  for 
formulation  of  corrections  to  atmospheric  density  models. 

2.2  Methodology 

The  results  presented  in  this  paper  were  generated  by  processing  the  positions  and  velocities  from  the 
CHAMP  and  GRACE  POE  data  as  observations  in  an  optimal  orbit  determination  process  in  order  to 
estimate  density  and  ballistic  coefficient.  The  POE  data  are  available  as  rapid  science  orbits  (RSO)  and 
may  be  downloaded  from  the  website  http://isdc.gfz-potsdam.de.  Konig  et  al.  (2002,  2005,  2006)  and 


Michalek  et  al.  (2003)  discuss  the  processing  and  accuracy  of  the  RSOs.  The  published  accuracy  of  the 
RSOs  compared  to  SLR  data  is  5-10  cm. 

Atmospheric  density  is  estimated  as  a  correction  to  a  baseline  atmospheric  density  model  as  part  of  an 
orbit  determination  scheme  using  Orbit  Determination  Tool  Kit  (ODTK).  The  POE  data  are  input  as 
measurements  into  a  sequential  measurement  processing  and  filtering  scheme.  A  smoother  is  then 
applied  to  the  filtered  solution  to  account  for  all  available  solution  data  to  increase  the  accuracy  over  the 
whole  solution  span.  The  filter  and  smoother  combination  estimates  the  time  variable  density  and 
ballistic  coefficient  including  realistic  covariance  matrices  established  by  the  physics  associated  with  the 
problem.  A  90x90  GRACE  Gravity  Model  2  (GGM02C)  (Tapley  et  al.,  2005),  solar  radiation  pressure, 
Earth  infrared  and  albedo  radiation  pressure,  luni-solar  point  masses,  general  relativity,  and  solid  Earth 
and  ocean  tides  are  force  models  in  addition  to  drag  included  in  the  orbit  determination  process.  SRP  was 
modeled  assuming  CHAMP  and  GRACE  were  spheres  with  area  of  6.5  m2,  which  represents  an  averaged 
area.  Earth  infrared  and  albedo  radiation  pressure  models  assumed  a  constant  area  of  20  m2,  where  the 
difference  from  the  SRP  area  accounts  for  the  larger  area  projected  toward  Earth  for  CHAMP  and 
GRACE.  Perfect  absorption  was  assumed  as  a  baseline  for  the  radiation  pressure  modeling,  but  the 
coefficients  are  estimated.  SRP  is  about  two  orders  of  magnitude  lower  than  drag  for  the  CHAMP 
altitude  so  these  simplified  models  should  introduce  minimal  error  (Bruinsma  et  al.,  2004).  Process  noise 
was  only  included  for  the  gravity  model.  Since  covariance  information  is  not  available  for  GGM02C,  the 
process  noise  level  was  determined  using  Joint  Gravity  Model  2  (JGM2)  (Nerem  et  al.,  1994)  error  levels 
for  low  Earth  orbit  satellites.  This  is  the  default  setting  in  ODTK.  This  is  a  conservative  process  noise 
level  for  the  gravity  model  since  the  JGM2  errors  are  larger  than  the  GGM02C  errors.  The  POE 
observations  were  provided  every  30  s  and  given  a  noise  level  of  10  cm. 

Wright  (2003)  and  Wright  and  Woodburn  (2004)  outline  techniques  for  estimating  density  that  are 
available  in  ODTK  software  package,  which  was  used  for  this  research.  The  technique  permits  the  local 
atmospheric  density  to  be  estimated  in  real-time  in  conjunction  with  the  orbit  determination  process  and 
provides  a  significant  improvement  over  the  standard  technique  of  estimating  the  ballistic  coefficient 
(BC)  or  drag  coefficient  because  BC  estimates  absorb  the  errors  generated  by  the  model  for  atmospheric 
density  and  BC.  Additionally,  BC  estimates  will  often  include  geopotential  model  errors.  Tanygin  and 
Wright  (2004)  present  the  polynomial  spline  fit  technique  for  geomagnetic  indices  utilized  in  the  orbit 
determination  scheme. 

The  estimated  atmospheric  density  is  obtained  as  a  correction  to  an  atmospheric  model.  The  models 
available  in  ODTK  are  Jacchia  1971  (Jacchia,  1971),  Jacchia-Roberts  (Roberts,  1971),  CIRA  1972 
(COSPAR  Working  Group  IV,  1972),  MSISE-1990  (Hedin,  1991),  and  NRLMSISE-2000  (Picone  et  al., 
2002).  The  Jacchia-Roberts  and  CIRA  1972  density  models  were  derived  from  the  Jacchia  1971  model 
so  the  results  from  these  three  models  should  be  similar  and  these  will  be  classified  as  the  Jacchia  models. 
Once  a  baseline  atmospheric  density  model  is  selected,  two  different  types  of  corrections  are  applied  to 
the  model  as  part  of  the  orbit  determination  process.  The  first  is  a  baseline  correction  derived  from 
historical  F10.7  and  ap  measurements  obtained  over  several  solar  cycles.  The  baseline  density  correction 
propagates  from  perigee  height  through  the  use  of  an  exponential  Gauss-Markov  sequence.  Error  in  the 
atmospheric  density  at  perigee  is  related  to  the  current  point  in  the  orbit  through  a  transformation  internal 
to  ODTK. 

The  second  type  of  atmospheric  density  correction  utilizes  the  observations  and  current  conditions  to 
yield  a  dynamic  correction.  The  use  of  a  sequential  process  allows  a  correction  to  be  estimated  at  each 
time  step  in  the  filter  as  opposed  to  a  single  correction  applied  to  the  complete  time  span  of  data  or  many 
corrections  increasing  the  size  of  the  state  for  a  batch  least  squares  process. 


Exponentially  correlated  Gauss-Markov  processes  describing  the  modeling  errors  also  exist  for  dynamic 
corrections  as  observed  for  the  baseline  corrections.  The  user  of  the  optimal  orbit  determination  process 


may  specify  the  values  of  the  density  and  ballistic  coefficient  exponential  Gauss-Markov  process  half- 
lives,  which  determines  the  effect  of  past  data  on  the  individual  density  correction  at  each  time  step. 
While  the  ballistic  coefficient  is  estimated  as  part  of  the  filter/smoother  process,  the  ballistic  coefficient 
was  initialized  using  yearly  averages  of  0.00444  m2/kg  for  2002-2003  and  0.00436  m2/kg  for  2004-2005 
(Bowman  et  al.,  2008).  Values  for  the  CHAMP  satellite’s  nominal  ballistic  coefficient  that  were  not 
included  in  these  years  were  estimated  for  years  both  preceding  and  following  these  ranges  by  taking  into 
account  the  changing  mass  of  the  satellite  (Hiatt,  2009).  The  nominal  ballistic  coefficient  for  GRACE  is 
defined  as  0.00687  m2/kg  in  the  ODTK  orbit  determination  scheme  (Bowman  et  al.,  2008). 

The  derived  densities  obtained  from  the  POE  data  are  compared  against  those  derived  from  the  CHAMP 
accelerometer  as  calculated  by  Sean  Bruinsma  of  CNES.  Bruinsma’s  density  values  are  averaged  to  10- 
second  intervals.  POE  derived  density  solutions  were  developed  for  multiple  time  periods  with  varying 
levels  of  solar  and  geomagnetic  activity.  The  authors  have  used  geomagnetic  and  solar  activity  bins 
defined  based  on  Picone  et  al.  (2002)  to  establish  low  solar  activity  as  defined  by  F10.7  <  75,  moderate 
solar  activity  for  75  <  F10.7  <  150,  elevated  solar  activity  for  150  <  F10.7  <  190,  and  high  solar  activity 
by  F10.7  >  190.  Additionally,  quiet  geomagnetic  periods  are  defined  for  Ap  <  10,  moderate  geomagnetic 
activity  for  10  <  Ap  <  50,  and  active  geomagnetic  periods  by  Ap  >  50.  These  definitions  are  applied  to 
separate  the  POE  derived  density  solutions  into  the  appropriate  categories.  This  cataloging  of  solutions 
facilitates  the  analysis  of  the  POE  derived  density  solutions  as  a  function  of  solar  and  geomagnetic 
activity. 

2.3  Initial  CHAMP  Results  and  Calibration 

McLaughlin  and  Bieber  (2007)  presented  the  initial  results  of  estimating  density  using  CHAMP  POE. 

The  paper  found  densities  estimated  using  different  baseline  density  models  and  densities  estimated  in 
overlap  regions  to  be  consistent  within  10  percent.  McLaughlin  et  al.  (2008a)  made  the  first  comparisons 
of  CHAMP  POE  derived  atmospheric  densities  to  those  derived  from  the  accelerometer  data  and  found 
errors  bounded  within  about  10  percent  for  days  with  moderate  solar  and  geomagnetic  activity.  Follow- 
on  work  focused  on  calibrating  the  POE  derived  densities  using  the  accelerometer  derived  densities. 
McLaughlin  et  al.  (2008b,  201  lb),  Hiatt  et  al.  (2009),  and  Hiatt  (2009)  presented  initial  calibration  results 
for  CHAMP  and  showed  that  density  and  ballistic  coefficient  half-lives  of  1 800  and  1 80,000  minutes  did 
not  perform  well  in  estimating  density.  Therefore,  further  work  focused  on  half-lives  of  1 .8,  18,  or  180 
minutes.  These  works  also  showed  that  POE  derived  corrections  to  Jacchia  based  density  models 
matched  the  accelerometer  derived  densities  better  than  POE  derived  corrections  to  MSIS  based  density 
models  did.  Lechtenberg  (2010)  included  a  more  comprehensive  examination  of  CHAMP  POE  derived 
densities  and  examined  days  from  2001  to  2007  picked  so  that  the  distribution  of  solar  and  geomagnetic 
activity  level  bins  matched  those  over  all  days  from  the  launch  of  CHAMP  through  2007  and  so  that  days 
from  different  times  in  the  year  were  adequately  represented.  The  results  of  this  study  are  summarized  in 
Table  1  including  the  best  density  half-life,  ballistic  coefficient  half-life,  and  baseline  density  model 
combination.  The  best  combination  was  determined  based  on  both  cross  correlation  and  RMS.  The 
different  Jacchia  based  density  models  had  nearly  identical  results,  especially  for  cross  correlation. 

These  results  show  the  very  high  correlation  between  CHAMP  POE  and  accelerometer  derived  densities 
and  relatively  low  RMS  differences.  The  consistently  best  combination  included  a  ballistic  coefficient 
half-life  of  1.8  minutes  and  a  CIRA  1972  baseline  density  model.  In  addition,  180  minutes  was  usually 
the  best  half-life  for  density.  The  differences  between  density  found  with  1 8  and  1 80  minute  half-lives 
for  density  was  usually  negligibly  different. 


Table  1 :  Best  Combinations  for  CHAMP  POE-Derived  Density  Correlation  and  RMS  with 


Accelerometer-Derived  Density. 


Bin 

Density  Half- 
Life  (min) 

BC  Half-Life 
(min) 

Atmospheric 

Model 

Cross 

Correlation 

RMS 

(10  12  kg/m3) 

Overall 

180 

1.8 

CIRA  1972 

0.910 

0.570 

Quiet  Geomagnetic 

180 

1.8 

CIRA  1972 

0.955 

0.348 

Moderate  Geomagnetic 

180 

1.8 

CIRA  1972 

0.919 

0.464 

Active  Geomagnetic 

18 

1.8 

CIRA  1972 

0.856 

0.960 

Low  Solar 

180 

1.8 

CIRA  1972 

0.931 

0.296 

Moderate  Solar 

180 

1.8 

CIRA  1972 

0.893 

0.559 

Elevated  Solar 

18 

1.8 

CIRA  1972 

0.947 

0.581 

Active  Solar 

180 

1.8 

CIRA  1972 

0.906 

0.759 

2.4  CHAMP  and  GRACE  Calibration 

Fattig  et  al.  (2010),  Fattig  (201 1),  and  McLaughlin  et  al.  (2013)  presented  the  first  extensive  calibration 
that  examined  the  same  days  for  CHAMP  and  GRACE  to  see  if  the  best  parameters  were  consistent 
between  CHAMP  and  GRACE.  These  studies  included  days  from  the  launch  of  GRACE  in  2004  through 
2009.  No  elevated  or  high  solar  activity  days  occurred  during  this  time  period.  Table  2  shows  a  summary 
of  the  results  for  CHAMP  and  Table  3  shows  a  summary  of  the  GRACE  results. 

Tables  2  and  3  show  that  the  results  are  fairly  consistent  between  the  CHAMP  and  GRACE  orbits  and 
across  solar  and  geomagnetic  activity  levels.  Again  the  overall  best  combination  was  with  corrections  to 
a  Jacchia-based  model.  Ballistic  coefficient  half-life  was  best  at  1.8  minutes  and  the  density  half-life  was 
best  at  18  or  180  minutes  with  the  difference  between  the  densities  estimated  with  the  two  density  half- 
lives  being  small. 

Table  2:  Best  Combinations  for  CHAMP  POE-Derived  Density  Correlation  and  RMS  with 


Accelerometer-Derived  Density. 


Bin 

Density  Half- 
Life  (min) 

BC  Half- 
Life  (min) 

Atmospheric 

Model 

Cross 

Correlation 

RMS 

(10' 12  kg/m3) 

Overall 

18 

1.8 

CIRA  1972 

0.883 

0.623 

Quiet  Geomagnetic 

180 

1.8 

CIRA  1972 

0.936 

0.272 

Moderate  Geomagnetic 

180 

1.8 

CIRA  1972 

0.945 

0.418 

Active  Geomagnetic 

18 

1.8 

CIRA  1972 

0.818 

0.961 

Low  Solar 

180 

1.8 

CIRA  1972 

0.928 

0.264 

Moderate  Solar 

18 

1.8 

CIRA  1972 

0.868 

0.761 

Table  3:  Best  Combinations  for  GRACE  POE-Derived  Density  Correlation  and  RMS  with 

Accelerometer-Derived  Density. 


Bin 

Density 

Half-Life 

(min) 

BC  Half- 
Life  (min) 

Atmospheric 

Model 

Cross 

Correlation 

RMS 

(1012  kg/m3) 

Overall 

180 

1.8 

Jacchia-Roberts 

0.849 

0.123 

Quiet  Geomagnetic 

180 

1.8 

Jacchia  1971 

0.802 

0.037 

Moderate  Geomagnetic 

180 

1.8 

Jacchia-Roberts 

0.912 

0.076 

Active  Geomagnetic 

180 

1.8 

Jacchia-Roberts 

0.845 

0.207 

Low  Solar 

180 

1.8 

Jacchia  1971 

0.766 

0.030 

Moderate  Solar 

180 

1.8 

Jacchia-Roberts 

0.883 

0.160 

McLaughlin  et  al.  (2009)  and  Mehta  and  McLaughlin  (201 1)  examined  the  simultaneous  ballistic 
coefficient  and  density  estimation.  The  estimated  density  magnitude  is  sensitive  to  the  nominal  ballistic 
coefficient  selected  so  a  bias  in  nominal  ballistic  coefficient  does  lead  to  a  bias  in  the  estimated  density. 
Therefore,  ballistic  coefficient  and  density  are  not  completely  separable  without  a  good  knowledge  of 
what  the  average  ballistic  coefficient  should  be.  The  estimated  ballistic  coefficients  show  a  3-5%  percent 
variation  in  magnitude  with  a  cycle  of  roughly  twice  per  orbit.  For  CHAMP,  the  correlation  between  the 
density  error,  defined  as  the  difference  between  accelerometer  and  POE  derived  densities,  and  the  ballistic 
coefficient  was  around  -0.3  for  ballistic  coefficient  half-life  of  1.8  min  and  density  half-life  of  18  min  and 
near  zero  if  density  half-life  was  increased  to  1 80  min.  For  GRACE  the  correlation  for  the  density  half- 
life  of  1 80  min  was  0. 1  with  other  results  about  the  same  as  for  CHAMP.  If  ballistic  coefficient  error  was 
being  absorbed  by  the  density  estimation,  one  would  expect  a  high  negative  correlation  between  these  two 
quantities.  However,  some  ballistic  coefficient  error  is  almost  certainly  being  absorbed  in  the  density 
estimates,  but  the  effects  on  the  density  estimation  are  small  based  on  comparisons  with  the 
accelerometer.  These  papers  also  examined  the  effects  of  fixing  ballistic  coefficient  at  the  nominal  value 
and  showed  that  the  difference  in  density  is  small  between  estimating  and  not  estimating  the  ballistic 
coefficient. 

Based  on  the  density  and  ballistic  coefficient  studies,  nominal  settings  were  used  for  further  analysis.  The 
nominal  settings  were  CIRA  1972  as  the  baseline  density  model,  ballistic  coefficient  half-life  of  1.8 
minutes,  and  density  half-life  of  180  minutes.  An  additional  benefit  of  180  minute  half-life  compared  to 
the  1 8  minute  half-life  is  that  it  provides  an  additional  order  of  magnitude  difference  between  the  ballistic 
coefficient  and  density  half-lives.  Large  separation  between  density  and  ballistic  coefficient  half-lives 
was  a  condition  for  the  two  factors  to  be  simultaneously  observable  (Wright,  2003;  Wright  and 
Woodburn,  2004). 

2.5  Continuous  Density  Data  Sets  and  Accuracy  Comparisons 

Mysore  Krishna  and  McLaughlin  (2011)  and  Mysore  Krishna  (2012)  presented  the  technique  used  to 
combine  14  hour  density  solutions  into  a  continuous  density  data  set  over  the  life  of  CHAMP  and  from 
launch  through  2011  for  GRACE.  The  data  are  available  as  one  week  data  sets  to  keep  the  individual  file 
sizes  manageable.  The  method  used  to  merge  the  data  is  the  linear  weighted  blending  technique  (Kim, 
2008;  Kim  et  al.,  2008;  Kim  et  al.,  2009).  The  technique  applies  weights  to  overlapping  areas  of  the 
fourteen  day  density  solution  based  on  how  far  the  point  is  from  the  end  or  beginning  of  the  overlapping 
solution.  The  14  hour  solutions  have  two  hour  overlaps  so  at  the  beginning  of  the  two  hour  overlap  the 
combined  density  will  match  the  density  of  the  earlier  of  the  two  fourteen  hour  solutions;  at  the  midpoint 
of  the  overlap  region  the  combined  density  will  be  a  straight  average  of  the  two  densities;  and  at  the  end 
of  the  overlap  region  the  combined  density  will  match  the  density  of  the  latter  of  the  two  fourteen  hour 
solutions.  The  primary  value  of  this  technique  is  to  provide  a  smooth  solution,  but  it  also  gives  more 
weight  to  points  farthest  from  the  endpoint  of  a  given  solution.  This  is  valuable  since  for  an  orbit  solution 
one  would  typically  expect  the  points  at  the  end  a  given  solution  to  be  of  lower  accuracy  than  points 
closer  to  the  middle.  Effectively,  this  technique  uses  information  from  the  other  overlapping  solution  to 
provide  more  information  about  what  is  happening  in  the  overlap  regions. 

Cross  correlation  and  RMS  of  the  HASDM,  POE  derived,  NRLMSISE-00,  and  Jacchia  71  densities  with 
accelerometer  derived  densities  were  calculated  and  compared  for  both  CHAMP  and  GRACE.  Table  4 
contains  the  cross  correlation  coefficients  of  the  POE  derived  density,  HASDM  densities,  and  of  the 
Jacchia  71  model  compared  to  accelerometer  derived  density  for  the  various  levels  of  solar  and 
geomagnetic  activity  as  well  as  for  the  entire  mission  life.  The  RMS  comparisons  are  contained  in  Table 
5.  The  tables  show  the  POE  derived  density  matches  the  accelerometer  density  significantly  better  than 
Jacchia  71  and  NRLMSISE-00  both  overall  and  for  every  level  of  solar  and  geomagnetic  activity.  The 
POE  derived  density  is  also  moderately  better  than  the  HASDM  density  for  all  conditions.  CHAMP  is 


included  as  a  satellite  in  the  HASDM  solutions  so  the  HASDM  densities  for  CHAMP  are  likely  better 
than  for  a  satellite  not  included  in  the  HASDM  solutions. 


Table  4:  Cross  correlation  coefficients  for  POE  derived  density,  HASDM,  NRLMSISE-00,  and  Jacchia 

7 1  with  accelerometer  derived  density  along  the  CHAMP  orbit. 

Bin 

Cross  Correlation 

POE  derived  density 

HASDM 

NRLMSISE-00 

Jacchia-7 1 

Overall 

0.934 

0.924 

0.888 

0.886 

Low  Solar 

0.926 

0.910 

0.880 

0.884 

Moderate  Solar 

0.935 

0.925 

0.881 

0.884 

Elevated  Solar 

0.938 

0.936 

0.907 

0.895 

High  Solar 

0.948 

0.942 

0.903 

0.895 

Quiet  Geomagnetic 

0.935 

0.923 

0.891 

0.896 

Moderate  Geomagnetic 

0.932 

0.924 

0.885 

0.874 

Active  Geomagnetic 

0.950 

0.941 

0.871 

0.831 

Table  5:  RMS  for  POE  derived  density,  HASDM,  NRLMSISE-00,  and  Jacchia  71  with  accelerometer 

derived  density  along  the  CHAMP  orbit. 

Bin 

Root  Mean  Square  (10 

12  kg/mJ) 

POE  derived  density 

HASDM 

NRLMSISE-00 

Jacchia-7 1 

Overall 

0.383 

0.400 

0.701 

Low  Solar 

0.322 

0.346 

0.849 

Moderate  Solar 

0.354 

0.372 

0.643 

Elevated  Solar 

0.526 

0.531 

0.719 

High  Solar 

0.573 

0.576 

0.783 

1.434 

Quiet  Geomagnetic 

0.346 

0.361 

0.745 

0.759 

Moderate 

Geomagnetic 

0.423 

0.441 

0.658 

0.895 

Active  Geomagnetic 

0.925 

0.986 

0.551 

3.231 

Table  6  shows  the  cross  correlation  and  Table  7  shows  the  RMS  results  for  GRACE.  No  elevated  or  high 
solar  activity  days  or  active  geomagnetic  weeks  occurred  during  the  available  GRACE  POE  data  (stalling 
in  2004).  The  trends  are  very  similar  to  those  observed  in  the  CHAMP  data.  Once  again  the  POE  derived 
density  is  clearly  superior  to  the  Jacchia  71  and  NRLMSISE-00  models  and  slightly  better  than  HASDM 
in  matching  the  accelerometer  derived  density.  These  results  show  the  promise  of  POE  derived  density 
from  many  satellites  to  improve  the  overall  knowledge  of  density  in  the  thermosphere. 


Table  6:  Cross  correlation  coefficients  for  POE  derived  density,  HASDM,  NRLMSISE-00,  and  Jacchia 
71  with  accelerometer  derived  density  along  the  GRACE  orbit. 


Bin 

Cross  Correlation 

POE  derived  density 

HASDM 

NRLMSISE-00 

Jacchia-7 1 

Overall 

0.885 

0.873 

0.844 

0.839 

Low  Solar 

0.855 

0.822 

ansa 

Moderate  Solar 

0.912 

0.863 

— 

Quiet  Geomagnetic 

0.883 

0.869 

0.852 

■SB 

Moderate  Geomagnetic 

0.891 

0.881 

0.827 

wm.m 

Table  7:  RMS  for  POE  derived  density,  HASDM,  NRLMSISE-00,  and  Jacchia  71  with  accelerometer 

derived  density  along  the  GRACE  orbit. 


Bin 

Root  Mean  Square  (10 

12  kg/m1) 

POE  derived  density 

HASDM 

NRLMSISE-00 

Jacchia-71 

For  all  Bins 

0.044 

0.047 

0.089 

0.111 

Low  Solar 

0.031 

0.031 

0.079 

0.100 

Moderate  Solar 

0.056 

0.060 

0.098 

0.122 

Quiet  Geomagnetic 

0.036 

0.038 

0.080 

0.096 

Moderate  Geomagnetic 

0.063 

0.068 

0.109 

0.147 

Figure  1  shows  the  cross  correlation  of  CHAMP  POE  derived  and  HASDM  densities  with  accelerometer 
derived  densities  and  the  angle  between  the  CHAMP  orbit  plane  and  the  Sun  vector,  P,  from  2001  to 
2008.  Figure  2  shows  the  same  for  GRACE.  Both  figures  show  a  periodic  variation  in  correlation  for 
both  POE  and  HASDM  densities  where  the  lower  correlation  occurs  where  P  is  near  zero.  This  occurs 
when  the  orbit  plane  is  near  the  terminator.  The  reason  for  the  lower  correlation  is  that  the  normal 
daylight/darkness  variation  that  dominates  the  signal  is  much  lower  in  amplitude  when  the  orbit  plane  is 
near  the  terminator.  When  this  happens  the  high  frequency  density  variations  that  are  only  observed  by 
the  accelerometer  are  of  the  same  order  as  the  daylight/darkness  variations,  which  brings  the  correlations 
down.  This  is  discussed  more  fully  in  McLaughlin  et  al.  (201  la). 


Figure  1 :  Beta  angle  and  weekly  cross  correlation  between  POE  and  HASDM  densities  with 

accelerometer  densities  for  CHAMP. 


Figure  2:  Beta  angle  and  weekly  cross  correlation  between  POE  and  HASDM  densities  with 

accelerometer  densities  for  GRACE. 


2.6  Comparison  of  Density  Estimation  from  Multiple  Satellites 

Interesting  times  to  compare  the  density  found  from  the  CHAMP  and  GRACE  satellites  are  when  the 
satellites’  orbits  are  nearly  coplanar.  One  such  time  occurred  during  April  of  2005.  Of  particular  interest 
is  when  GRACE  is  flying  nearly  directly  above  CHAMP  because  the  density  differences  can  be  attributed 
solely  to  altitude.  For  the  April  2005  coplanar  period  this  condition  was  met  for  a  few  hours  about  every 
two  days  including  on  April  3  and  April  5.  Solar  activity  was  at  moderate  levels  on  both  days. 
Geomagnetic  activity  was  low  on  April  3  and  right  on  the  moderate/active  border  on  April  5.  However,  a 
geomagnetic  storm  occurred  on  late  April  4,  extending  into  early  April  5  where  ap  peaked  at  132  early  on 
April  5.  This  geomagnetic  activity  is  apparent  in  the  data  for  April  5.  The  analysis  used  CIRA  1972  as  the 
baseline  density  model,  BC  half-life  of  1.8  min,  and  density  half-life  of  180  min  for  both  CHAMP  and 
GRACE  POE  densities. 

Figure  3  shows  the  densities  for  CHAMP  and  GRACE- A  for  April  3,  2005  found  from  the  accelerometer, 
from  HASDM,  from  Jacchia  1971,  from  precision  orbit  ephemerides  (POE),  and  from  NRLMSISE-00. 
Because  of  the  higher  altitude,  GRACE  densities  are  an  order  of  magnitude  lower  than  the  CHAMP 
densities,  but  the  overall  structure  is  similar  between  the  satellites  for  each  individual  data  set.  Interesting 
features  observable  in  the  accelerometer  data  are  nocturnal  peaks  in  density  that  occur  for  both  satellites 
between  the  major  diurnal  peaks  that  occur  during  the  sunlit  portion  of  the  orbit.  These  nocturnal  features 
are  not  observable  in  most  of  the  other  data  sources  except  for  some  very  slight  humps  or  elevated  density 
troughs  during  the  times  of  the  nocturnal  peaks.  However,  NRLMSISE-00  does  model  these  peaks,  but 
for  GRACE  often  has  a  secondary  peak  that  does  not  appear  in  the  accelerometer  data.  The  nocturnal 
peaks  occur  over  the  equator  and  correspond  to  the  Midnight  Density  Maximum  (MDM)  discussed  in 
Arduini  et  al.  and  Akmaev  et  al.  Both  the  POE  and  HASDM  densities  match  the  accelerometer  very  well 
for  the  CHAMP  satellite,  but  the  Jacchia  1971  and  NRLMSISE-00  models  tend  to  overestimate  peak 
density.  For  GRACE  all  other  data  sources  overestimate  the  sunlit  peaks  seen  in  the  accelerometer  data. 
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Figure  3.  Densities  Measured  and  Estimated  for  the  CHAMP  and  GRACE- A  Satellites  on  April 

3,  2005. 


Figure  4.  Densities  Measured  and  Estimated  for  the  CHAMP  and  GRACE- A  Satellites  on  April  5, 

2005. 


Figure  4  shows  data  similar  to  Figure  3  for  the  next  time  period  when  GRACE  is  nearly  directly  above 
CHAMP  on  April  5,  2005.  The  density  magnitudes  are  approximately  double  what  they  were  on  April  3. 
This  increase  is  caused  by  the  geomagnetic  storm  that  began  on  late  April  4  as  discussed  earlier.  The  early 
part  of  this  time  period  shows  an  overestimation  of  atmospheric  density  by  Jacchia  1971  caused  by  the 
elevated  geomagnetic  indices  earlier  in  the  day.  The  densities  for  both  CHAMP  and  GRACE  gradually 
decline  to  lower  values  over  the  course  of  this  time  period.  Again,  both  CHAMP  and  GRACE 
accelerometer  densities  observe  MDM  not  well  observed  in  the  other  data  sources  except  NRLMSISE-00. 
HASDM  and  the  POE  densities  both  match  the  accelerometer  densities  very  well  for  CHAMP,  but  not  as 
well  for  GRACE.  NRLMSISE-00  does  quite  well  for  this  time  period  except  for  the  extra  peaks 
appearing  for  GRACE. 

This  subsection  also  examines  densities  estimated  from  CHAMP,  GRACE,  and  TerraSAR-X  for  the 
period  of  September  21-30,  2007.  During  this  period,  CHAMP  had  an  altitude  of  about  360  km,  GRACE 
had  an  altitude  of  about  473  km,  and  TerraSAR-X  had  an  altitude  of  about  528  km.  Simultaneous  density 
estimation  for  the  three  satellites  shows  atmospheric  activity  at  multiple  altitudes  in  the  atmosphere. 

CIRA  1972  was  used  as  the  baseline  density  model  with  180  min  for  the  density  half-life  and  1.8  min  for 
the  ballistic  coefficient  half-life  for  all  satellites. 


Figure  5.  Estimated  and  Measured  Densities  for  CHAMP,  GRACE,  and  TerraSAR-X,  September  26- 

27,  2007. 

Figure  5  shows  the  POE,  Jacchia  1971,  and  NRLMSISE-00  densities  for  CHAMP,  GRACE,  and 
TerraSAR-X;  and  the  accelerometer  densities  for  CHAMP  and  GRACE  for  a  14  hour  solution  from 
September  26  to  September  27,  2007.  Solar  activity  level  was  low  on  both  days  and  geomagnetic  activity 
was  low  on  the  26th  and  moderate  on  the  27th.  CHAMP  densities  are  more  than  an  order  of  magnitude 
greater  than  GRACE  densities  and  TerraSAR-X  densities  are  about  half  of  GRACE  densities.  The  MDM 
seen  in  the  previous  section  are  apparent  in  the  CHAMP  accelerometer  data,  but  do  not  appear  in  the  POE 
density.  These  peaks  are  also  no  longer  readily  observable  in  the  GRACE  accelerometer  derived 
densities.  However,  GRACE  is  not  flying  nearly  directly  above  CHAMP  as  in  the  previous  section. 


Aside  from  these  peaks  and  some  other  short  period  variations,  the  POE  density  matches  the 
accelerometer  density  very  well.  NRLMSISE-00  and  Jacchia  1971  significantly  overestimated  the 
density  for  both  CHAMP  and  GRACE  compared  to  the  accelerometer  and  POE  density.  For  TerraSAR- 
X  the  POE  density  is  very  near  Jacchia  1971 .  This  could  be  caused  by  either  Jacchia  1971  better  modeling 
the  density  at  this  higher  altitude  or  the  POE  technique  having  increased  difficulty  observing  density 
variations  at  higher  satellite  altitudes  and  lower  density  magnitudes.  Although  not  shown  in  the  figure, 
the  POE  density  derived  as  corrections  to  NRLMSISE-00  for  TerraSAR-X  was  close  to  the  POE  density 
derived  from  CIRA  1972  that  is  shown  in  the  figure. 

Figure  6  shows  data  similar  to  Figure  5  for  14  hours  from  September  29  to  30,  2007.  Solar  activity  level 
was  low  on  both  days  and  geomagnetic  activity  was  moderate  on  the  29th  and  30th.  The  trends  are 
similar,  but  the  density  magnitudes  for  all  satellites  are  approximately  double  the  values  on  September  26 
and  27.  The  MDM  are  again  seen  in  the  CHAMP  accelerometer  densities,  but  are  not  well  defined  in 
GRACE.  Jacchia  1971  again  overestimates  density  for  CHAMP  and  GRACE,  but  POE  derived  densities 
match  the  accelerometer  very  well.  NRLMSISE-00  does  very  well  for  CHAMP,  but  overestimates  the 
densities  for  GRACE.  For  TerraSAR-X  the  difference  between  POE  and  Jacchia  1971  densities  is  again 
much  less  distinct  than  for  CHAMP  and  GRACE,  but  the  values  have  greater  separation  than  on 
September  26  and  27.  This  difference  could  be  caused  by  higher  overall  density  on  the  later  days. 
NRLMSISE-00  and  the  POE  densities  match  well  for  TerraSAR-X  for  the  last  half  of  this  time  span. 


Figure  6.  Estimated  and  Measured  Densities  for  CHAMP,  GRACE,  and  TerraSAR-X,  September 

29-30,  2007. 

2.7  Effects  of  Orbit  Ephemeris  Error  and  Limited  Data  on  Density  Estimation 

The  previous  work  has  focused  on  satellites  with  continuous  data  of  high  accuracy,  which  is  only 
available  from  a  limited  number  of  satellites  with  high  quality  GPS  receivers.  The  effects  of  limited  data 
and  lower  quality  data  on  density  estimation  are  examined  to  see  the  suitability  of  using  a  larger  set  of 


satellites  to  estimate  density.  Discontinuous  data  sets  are  created  by  decimating  the  original  POE  data  for 
CHAMP  and  GRACE.  This  file  is  then  used  as  an  input  in  ODTK  to  obtain  the  POE  derived  density 
estimates. 


The  next  step  is  to  introduce  different  levels  of  ephemeris  noise  to  the  POE  data  and  estimate  density 
using  the  lower  accuracy  data.  The  effects  reduced  accuracy  ephemeris  data  and  its  effect  on  density 
estimation  was  initially  examined  by  Locke  (2012).  Noise  levels  of  0.1,  0.5,  1,  10,  and  100  m  were  added 
to  the  data  and  data  was  decimated  so  that  data  existed  for  15  and  8  minutes  per  orbit. 


Tables  8  and  9  show  the  cross  correlation  results  with  15  minutes  of  data  per  orbit  for  CHAMP  and 
GRACE,  respectively.  As  expected,  the  cross  correlation  goes  down  as  the  noise  level  is  increased. 
However,  the  overall  accuracy  remains  reasonably  high  even  up  to  100  m  noise  levels.  The  results  for 
RMS  and  for  8  minutes  per  orbit  show  similar  decreases  in  accuracy/precision  and  even  with  only  eight 
minutes  of  data  per  orbit  the  density  estimates  remain  usable  for  even  relatively  high  noise  levels. 


Table  8.  Cross  correlation  between  POE  derived  density  estimates  and  Bruinsma’s  accelerometer  derived 
_ densities  for  the  CHAMP  orbit  using  various  noise  level  data  available  15  minutes  per  orbit _ 


Cross  Correlation  for  various  noise  levels  (m) 


Bin 

Full 

Data 

0 

0.1 

0.5 

1.0 

10.0 

100.0 

Overall 

0.917 

0.913 

0.906 

0.901 

0.899 

0.893 

0.869 

Low  Solar 

0.900 

0.888 

0.877 

0.852 

0.846 

0.847 

0.820 

Moderate 

Solar 

0.922 

0.920 

0.914 

0.914 

0.914 

0.907 

0.887 

Elevated  Solar 

0.948 

0.946 

0.942 

0.945 

0.944 

0.938 

0.919 

High  Solar 

0.890 

0.880 

0.877 

0.876 

0.876 

0.857 

0.811 

Quiet 

Geomagnetic 

0.937 

0.934 

0.927 

0.921 

0.920 

0.916 

0.900 

Moderate 

Geomagnetic 

0.892 

0.888 

0.881 

0.877 

0.875 

0.867 

0.829 

Active 

Geomagnetic 

0.855 

0.841 

0.838 

0.838 

0.838 

0.821 

0.784 

Table  9.  Cross  correlation  between  POE  derived  density  estimates  and  Bruinsma’s  accelerometer  derived 
densities  for  the  GRACE  orbit  using  various  noise  level  data  available  15  minutes  per  orbit 


Cross  Correlation  for  various  noise  levels  (m) 


Bin 

Full 

Data 

0 

0.1 

0.5 

1.0 

10.0 

100.0 

Overall 

0.862 

0.825 

0.820 

0.813 

0.813 

0.815 

0.813 

Low  Solar 

0.708 

0.598 

0.601 

0.587 

0.588 

0.603 

0.620 

Moderate 

Solar 

0.901 

0.886 

0.874 

0.875 

0.875 

0.871 

0.854 

Quiet 

Geomagnetic 

0.850 

0.808 

0.807 

0.794 

0.797 

0.795 

0.802 

Moderate 

Geomagnetic 

0.861 

0.827 

0.817 

0.823 

0.818 

0.826 

0.837 

Active 

Geomagnetic 

0.831 

0.813 

0.808 

0.787 

0.784 

0.776 

0.689 

3  Effects  on  Orbit  Propagation 
3.1  Background  and  Methodology 

Some  previous  work  has  examined  the  effects  of  different  frequency  signals  on  orbit  propagation. 
Anderson  et  al.  (2009)  studied  the  effects  of  different  wavelength  density  perturbations  and  temporal 
resolution  on  orbit  propagation.  Schaeperkoetter  and  McLaughlin  (2010)  compared  orbit  propagations 
along  the  CEIAMP  orbit  using  accelerometer  derived  density  to  those  using  POE  derived  density  and 
EIASDM  density  for  several  generic  time  periods.  However,  neither  of  these  works  directly  examines 
how  specific  types  of  high  frequency  signals  seen  in  accelerometer  derived  densities  affect  orbit 
propagation.  The  effects  have  clear  implications  for  the  need  to  model  these  variations  in  future  density 
models  used  in  orbit  analysis. 

McLaughlin,  Locke,  and  Mysore  Krishna  (2012)  use  a  precision  orbit  initial  vector  and  propagate  the 
orbit  using  central  body  and  J2  gravitation  with  drag  using  density  from  accelerometer,  POE,  HASDM,  or 
Jacchia  71.  The  orbits  using  the  different  densities  are  compared  for  24  hour  propagation  using  both  root 
mean  square  and  maximum  difference.  This  method  does  not  account  for  how  the  difference  in 
propagated  position  would  give  different  model  results.  This  is  justified  by  assuming  the  density 
difference  would  be  small. 

A  position  and  velocity  vector  is  taken  from  the  POE  solutions  to  initialize  the  integration.  The  force 
models  include  Earth  central  body  and  oblateness  (J2)  geopotential,  and  drag.  Orbits  are  propagated  with 
density  from  the  accelerometer,  POE,  HASDM,  or  Jacchia  71  and  the  resulting  orbits  are  compared 
assuming  the  accelerometer  derived  density  gives  the  truth  orbit.  Earth  constants  are  from  the  WGS-84 
model.  The  inverse  ballistic  coefficients  are  defined  to  be  0.00436  m2/kg  for  CHAMP  and  0.00687 
m2/kg  for  GRACE. 

The  POE,  HASDM,  and  Jacchia  71  densities  are  normalized  to  the  same  mean  as  the  accelerometer 
densities  because  the  primary  concern  in  this  research  is  the  effect  of  temporal  variations  on  the  orbit 
propagation  and  not  on  the  bias  between  the  density  data  sets  or  models.  The  normalization  is  performed 


by  dividing  the  POE  derived,  HASDM,  or  Jacchia  7 1  densities  by  their  mean  and  multiplying  by  the 
mean  of  the  accelerometer  derived  densities.  The  means  used  are  those  for  each  particular  24  hour  period 
of  propagation.  This  method  of  normalization  was  used  assuming  the  bias  between  the  data  sets  is  caused 
by  errors  in  the  inverse  ballistic  coefficients  used  in  solving  for  densities. 

Accelerometer  derived  densities  are  given  in  ten  second  increments.  The  integrator  step  size  was  also  ten 
seconds  to  ensure  that  the  integration  included  all  variations  in  density  observed  by  the  accelerometer. 

3.2  CHAMP  Orbit  Propagations 

Normal  Day.  Figure  7  shows  the  densities  and  orbit  propagation  results  for  a  normal  day  for  CHAMP. 

The  densities  in  the  figure  are  the  actual  model  or  derived  densities  that  have  not  been  normalized.  The 
reason  for  plotting  the  densities  before  normalization  is  so  the  different  density  sources  can  be  more  easily 
observed  in  the  figure.  The  POE  derived  densities  clearly  provide  better  orbit  propagation  accuracy  than 
HASDM  or  Jacchia  7 1 .  The  worst  case  error  for  the  propagation  with  POE  derived  densities  is  only  13 
m,  which  would  not  be  significant  for  most  applications. 


Time  (hours  since  00:00  UTC  on  10/27/2005) 

Figure  7.  Density  and  orbit  propagation  errors  along  the  CHAMP  orbit  on  a  normal  day  (October 

27,2005). 

Traveling  Atmospheric  Disturbance.  Figure  8  includes  densities  and  orbit  propagation  errors  for  a  day 
that  includes  a  TAD,  which  is  an  increase  in  density  caused  by  geomagnetic  activity  that  begins  at  high 
latitudes  and  propagates  toward  the  equator  resulting  in  constructive  interference  between  the  waves  near 
the  equator.  The  TAD  on  April  19,  2002  occurs  between  1 100  and  1500  UTC  and  the  resulting  density 
enhancements  can  be  observed  in  the  figure.  In  addition,  in  the  orbit  propagation  plot  the  direction  of  the 
drift  caused  by  errors  in  the  Jacchia  71  model  changes  at  about  this  time.  However,  the  errors  in  orbit 
propagation  using  the  POE  derived  densities  remain  quite  small.  Clearly,  the  high  frequency  signals 
observed  by  the  accelerometer,  but  not  in  the  POE  derived  densities  do  not  have  a  major  effect  on  orbit 
propagation. 
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Figure  8.  Density  and  orbit  propagation  errors  along  the  CHAMP  orbit  on  a  day  with  a  traveling 
atmospheric  disturbance  (TAD)  (April  19,  2002). 

Polar  Cusp  Density  Enhancements.  The  accelerometer  on  CHAMP  observed  density  enhancements  that 
occurred  around  the  polar  cusp.  The  causes  of  these  increases  in  density  have  been  an  active  area  of 
investigation.  These  are  seen  as  spikes  in  density  that  occur  on  both  sides  of  the  magnetic  pole,  but  are 
not  always  observable  and  not  necessarily  tied  to  increased  geomagnetic  activity.  Figure  9  shows 
densities  and  orbit  errors  for  March  21,  2002,  which  was  a  day  on  which  polar  cusp  enhancements  were 
observed.  A  small  spike  in  the  accelerometer  density  can  be  seen  at  this  time,  but  no  discernible  effects 
are  observed  on  any  of  the  orbit  propagations.  The  maximum  error  for  the  orbit  propagation  with  the 
POE  derived  densities  only  reaches  8  m. 
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Figure  9.  Density  and  orbit  propagation  errors  along  the  CHAMP  orbit  on  a  day  with  a  polar  cusp 

density  enhancement  (March  21,  2002). 
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Figure  10.  Density  and  orbit  propagation  errors  along  the  CHAMP  orbit  on  a  day  when  the 
CHAMP  orbit  plane  is  near  terminator  (January  31,  2002). 


Orbit  Plane  near  the  Terminator.  When  the  orbit  plane  is  near  the  terminator,  the  high  frequency 
variations  seen  in  the  accelerometer  derived  densities  can  be  of  the  same  magnitude  as  the  variations  from 
the  satellite  traveling  from  daylight  into  darkness  and  back.  This  is  because  the  normal  variations  are  of 
much  lower  magnitude  when  the  satellite  never  travels  into  full  nighttime  darkness  or  full  local  noon 
daylight.  Figure  10  shows  how  small  the  variations  are,  but  in  this  particular  case  the  high  frequency 
variations  are  of  relatively  small  magnitude  still.  The  effects  on  orbit  propagation  are  very  small  as  can 
be  seen  in  the  lower  plot  of  Figure  10. 


GRACE 


Normal  Day.  Figure  1 1  shows  the  densities  and  orbit  errors  for  a  normal  day  for  GRACE.  Notice  that 
GRACE  is  at  a  higher  altitude  than  CHAMP  for  most  of  its  mission  life  and  the  densities  are  about  an 
order  of  magnitude  lower.  Consequently,  the  orbit  errors  caused  by  density  errors  are  also  considerably 
lower  for  GRACE  than  for  CHAMP.  The  maximum  orbit  propagation  errors  here  are  all  less  than  10  m. 
However,  the  reader  should  remember  that  these  are  errors  for  densities  that  have  been  normalized  to  the 
same  mean.  The  Jacchia  71  model  would  clearly  have  higher  errors  without  the  normalization  because  of 
the  bias  between  the  Jacchia  7 1  densities  and  the  other  densities  that  appears  in  the  upper  plot  of  Figure 
11. 
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Figure  1 1 .  Density  and  orbit  propagation  errors  along  the  GRACE  orbit  on  a  normal  day  (October 

1,2005). 
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Figure  12.  Density  and  orbit  propagation  errors  along  the  CHAMP  orbit  on  a  day  when  the 
CHAMP  orbit  plane  is  near  terminator  (November  2,  2005). 


Orbit  Plane  near  the  Terminator.  The  effects  of  flying  near  the  terminator  are  more  readily  seen  in 
GRACE  results  than  in  those  for  CHAMP  as  seen  in  Figure  12.  The  daylight/darkness  variations  are  all 
very  small  and  the  high  frequency  signals  in  the  accelerometer  data  are  much  more  apparent  than  for 
CHAMP.  However,  the  effects  of  these  unmodeled  high  frequency  variations  on  the  orbit  propagation 
remain  quite  small  with  errors  for  HASDM  and  POE  densities  barely  exceeding  1  m. 


5  Satellite  Drag  Coefficients  and  Ballistic  Coefficients 


5.1  Methodology 

The  Direct  Simulation  Monte  Carlo  (DSMC)  method  is  a  computational  technique  for  the  simulation  of 
dilute  gas  flows  at  the  molecular  level  and  is,  to  date,  the  basic  numerical  method  in  the  kinetic  theory  of 
gases  and  rarefied  gas  dynamics.  The  DSMC  method  uses  a  cell  and  particle  approach  to  track 
representative  molecules,  while  probabilistically  selecting  candidates  for  intermolecular  collisions.  Every 
simulated  molecule  represents  W  molecules  of  real  gas,  where  W  is  the  statistical  weight  of  a  simulated 
molecule.  W  typically  lies  in  the  range  of  1012  to  1025  real  molecules  per  simulated  molecule. 

The  computational  domain  is  divided  into  mean  free  path  sized  cells  and  molecules  in  each  cell  are 
tracked  independently.  The  positions,  velocities,  and  internal  energies  (in  the  case  of  polyatomic 
molecules)  of  simulated  molecules  are  stored  in  memory  and  are  modified  with  time  as  the  molecules  are 
concurrently  tracked  through  representative  collisions  and  boundary  interactions  within  the  computational 
domain.  The  movement  and  collisions  of  molecules  are  decoupled  based  on  the  dilute  gas  approximation. 
At  each  time-step,  collision  probabilities  are  calculated  and  collisions  are  carried  out  only  between 
molecules  in  the  same  cell.  The  size  of  the  time-step  is  comparable  to  the  mean  collision  time15. 

The  satellite  or  spacecraft  geometry  is  inserted  into  the  flow  field  as  a  surface  mesh.  The  mesh  format 
differs  from  one  program  to  another.  Molecules  are  inserted  into  the  flow  field  either  through  initial 
creation  at  the  first  time  step,  a  surface  flux  from  the  boundaries,  or  the  reservoir  method.  In  typical 
DSMC  simulations,  like  the  flow  over  a  satellite  or  a  vehicle  in  the  Earth’s  atmosphere,  the  computational 
domain  is  part  of  a  larger  flow  environment.  The  boundaries  of  the  computational  domain  are  therefore 
often  set  to  be  free-stream  boundaries  where  molecules  are  allowed  to  leave  and  enter  the  computational 
domain  and  the  number  of  simulated  molecules  varies  with  time. 


The  Three-Dimensional  Visual  Program  (DS3V)  was  used  to  model  the  interactions  of  simulated 
molecules  with  parameterized  surfaces  and  intermolecular  collision  dynamics.  In  spite  of  the 
computationally  demanding  nature  of  DS3V,  it  was  chosen  for  two  main  reasons:  (i)  it  is  freely  available 
on  the  internet  (www.gab.com.au),  and  (ii)  it  is  highly  reliable  and  widely  used.  The  only  computational 
parameter  specified  by  the  user  is  the  initial  number  of  megabytes  used  for  storage.  The  size  of  the  cells 
used  to  discretize  the  computational  domain  is  set  as  a  function  of  this  initial  number  of  megabytes 
defined  by  the  user.  The  program  sets  all  other  computational  variables  automatically.  However,  an 
optional  menu  is  available  should  the  user  choose  to  define  the  computational  parameters  such  as  cell  size 
and  timestep. 

DS3V  does  not  allow  the  user  to  explicitly  specify  the  total  energy  accommodation  coefficient.  The 
default  is  complete  energy  accommodation.  A  workaround  for  this  limitation  of  DS3V  was  initially 
proposed  by  Pilinski  et  al.  (201 1);  however,  most  of  their  work  was  limited  to  complete  accommodation. 
In  order  to  simulate  partial  energy  accommodation  (with  diffuse  reemission)  with  DS3V,  the  temperature 
of  the  spacecraft  surface  is  set  equal  to  the  kinetic  temperature  of  the  reemitted  molecules.  This  forces  an 
assumption  of  single  reflection  for  individual  molecules.  This  means  that  once  a  molecule  interacts  with 
the  surface  of  the  spacecraft,  the  molecule  is  assumed  to  travel  a  large  distance  before  colliding  with 
another  incoming  molecule  and  loses  any  chance  of  re -interacting  with  the  surface.  Single  reflections  are 
dominant  in  free  molecular  flow  for  simple  convex  geometries.  However,  caution  needs  to  be  taken  when 
dealing  with  concave  and  complex  geometries.  Pilinski  et  al.  (2011)  concluded  that  the  difference  in  drag 
coefficients  computed  for  concave  satellites  using  single  and  multiple  reflections  is  on  the  order  of  1%. 


5.2  Model  Development 

The  steps  involved  in  the  development  of  the  GRACE  drag  coefficient  model  are  presented  in  this 
section. 

Step  1 :  The  detailed  geometry  and  surface  properties  for  surface  temperature  calculations  were  obtained 
from  the  product  specification  document  developed  for  GRACE.  Equations  were  used  to  calculate  the 
worst-case  surface  temperatures  for  the  day  and  night  cycles  to  be  340  K  and  210  K  respectively.  The 
detailed  geometry  and  mesh  of  the  GRACE  satellite  developed  for  DSMC  computations  is  shown  in 
Figure  12.  A  highly  detailed  mesh  with  approximately  20,000  elements  captures  all  the  features  needed 
for  the  development  of  a  high  fidelity  model.  The  CAD  geometry  was  used  to  calculate  the  projection 
areas  at  different  attitude  orientations.  The  geometry  was  projected  onto  a  plane  perpendicular  to  the  flow 
and  the  outline  of  the  projection  was  used  to  estimate  the  areas.  Since  the  attitude  of  the  GRACE  satellites 
is  controlled,  the  attitude  bounds  were  set  at  -3°  and  3°  for  pitch  (O)  and  0°  and  3°  for  sideslip  (P). 
Symmetry  was  assumed  about  the  X-Z  plane.  Areas  were  estimated  at  1°  intervals  and  are  given  in  Table 
10.  The  CAD  geometry  model  and  estimated  areas  were  validated  using  the  front  and  side  profile  areas 
provided  in  the  GRACE  product  specification  manual. 
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Figure  13:  Detailed  geometry  and  mesh  with  orientation  of  the  GRACE  satellite  for  DSMC  computations. 


Table  10:  GRACE  cross-sectional  areas  as  a  function  of  attitude 


O  (right) 

P  (down) 

3° 

2° 

1° 

0° 

-1° 

-2° 

-3° 

0° 

1.311 

1.207 

1.104 

1.004 

1.093 

1.195 

1.299 

1° 

1.313 

1.210 

1.113 

1.049 

1.102 

1.198 

1.302 

2° 

1.315 

1.224 

1.156 

1.094 

1.145 

1.212 

1.303 

3° 

1.335 

1.266 

1.199 

1.139 

1.188 

1.254 

1.323 

Step  2:  DSMC  simulations  were  performed  for  all  attitude  orientations  at  a  large  sample  of  data  points 
along  the  orbit  of  GRACE.  Days  representing  different  solar  and  geomagnetic  activity  levels  were 
carefully  chosen  for  the  simulations.  Since  the  mission  life  of  GRACE  does  not  cover  elevated  or  high 
solar  activity  days,  simulated  solar  max  conditions  were  used.  The  atmospheric  properties  (atmospheric 
temperature,  total  number  density,  and  the  partial  pressures  of  the  atmospheric  species  constituents)  for 
the  DSMC  runs  were  calculated  using  the  NRLMSISE-00  model  inbuilt  into  the  Aerospace  Toolbox  in 
Matlab®.  The  semi-empirical  energy  accommodation  coefficient  model  developed  by  Pilinski  et  al.  (2010) 
was  used  with  all  its  assumptions  and  limitations.  The  energy  accommodation  coefficient  model  assumed 
a  diffuse  gas-surface  interaction  model. 


Step  3:  Since  DSMC  is  computationally  intensive,  CD  cannot  be  computed  real  time  and  a  semi -empirical 
model  is  needed  to  estimate  CD  at  any  point  in  the  orbit.  A  sensitivity  analysis  was  performed  on  the 
computed  CD  and  input  parameters  for  highest  correlations  to  develop  the  empirical  relations.  No  strong 
correlations  were  obtained  between  the  computed  total  CD  and  the  input  parameters.  However,  very  high 
correlations  were  seen  between  the  pressure  and  shear  components  of  CD  with  accommodation 
coefficient,  and  free-stream  atmospheric  temperature  respectively,  as  given  in  Table  11.  Based  on  the 
sensitivity  profiles  of  CD  with  the  different  input  parameter  obtained  in  Ref  20,  power  relations  were  used 
for  the  empirical  fits. 


Table  11:  Correlation  coefficients  of  pressure  and  shear  contributions  with  various  input  parameters. 


TX 

V„i 

XHe 

A0 

a 

Toted 

0.400 

-0.703 

-0.067 

-0.518 

-0.361 

Pressure 

-0.667 

-0.174 

0.698 

0.307 

-0.997 

Shear 

0.963 

-0.367 

-0.737 

-0.699 

0.705 

Step  4:  The  developed  model  relations  were  validated  using  additional  DSMC  simulations  that  were 
performed  at  random  points  in  time  along  the  orbit  of  GRACE.  The  error  in  CD  estimated  at  the  random 
points  using  the  model  relations  and  computed  using  DSMC  was  on  the  order  of  1%  with  a  mean  of  0.8%. 


5.3  Energy  Accommodation  Coefficient  Models 


The  semi-empirical  energy  accommodation  coefficient  (EC AM)  model  developed  by  Pilinski  et  al.  (2010) 
was  used  for  the  development  of  the  GRACE  CD  model.  The  limitation  of  the  energy  accommodation 
model  used  is  that  the  model  is  not  bounded  at  the  lower  end  by  the  energy  accommodation  coefficient  of 
the  substrate,  as  (clean  surface  with  no  adsorption  of  atomic  oxygen  to  the  surface),  originally  derived  by 
Goodman  (1966).  The  fraction  of  the  surface  covered  by  adsorbed  atomic  oxygen  becomes  smaller  as  the 
partial  pressure  of  atomic  oxygen  decreases.  Energy  accommodation  coefficient  of  a  substrate  material, 
as,  derived  in  Goodman  (1966)  is  given  as 


a,  = 


(1  +  ^)2 


Ks  is  the  substrate  coefficient,  and  p  is  the  mass  ratio  of  the  incoming  molecule  to  the  molecules  that 
compose  the  surface  lattice.  For  a  flat  plate  normal  to  the  flow,  Ks  is  2.4  when  the  distribution  of  the 
direction  of  the  incoming  molecules  is  random  (low  speed  ratio)  and  3.6  for  normal  incoming  flow  (high 
speed  ratio).  The  distribution  of  the  incoming  molecules  is  also  a  function  of  the  surface  geometry.  A 
sphere  has  a  distribution  of  incoming  molecules  that  is  random,  even  at  high  speed  ratios.  GRACE  has  a 
geometry  that  allows  normal  incidence  at  high  speed  ratios.  Since  the  attitude  of  GRACE  is  controlled 
and  the  computed  speed  ratios  lie  in  the  range  of  6  to  1 1  (relatively  high),  a  value  of  3.6  for  Ks  is  more 
suitable  for  calculating  the  accommodation  coefficient  of  the  substrate  for  GRACE. 


Modified  Semi-Empirical  Energy  Accommodation  Coefficient  Model 


The  modified  semi-empirical  energy  accommodation  coefficient  model  (EACM  (MOD))  is  a  modified 
version  of  the  original  model  with  Goodman’s  substrate  accommodation  implemented.  The 
implementation  of  Goodman’s  formula  is  achieved  in  a  very  simple  manner.  The  energy  accommodation 
coefficient  in  the  original  model  is  given  as 


1  +  K-n0-T„ 

Where  K  is  the  Langmuir  fitting  parameter  and  n0  is  the  number  density  of  atomic  oxygen.  The  energy 
accommodation  coefficient  with  Goodman’s  formula  implemented  is  given  as 


a  = 


a^  +  K  ■  rip  ■  T, 

1+  K  ■  n0  -Tx 


Semi-Empirical  Satellite  Accommodation  Model 

The  Semi-Empirical  Satellite  Accommodation  Model  (SESAM)  is  a  version  of  the  original  energy 
accommodation  coefficient  model  modified  by  Pilinski  et  al.  (2010).  The  modification  is  the 
implementation  of  Goodman’s  formula  to  the  original  model,  EC  AM.  Energy  accommodation 
coefficients  as  given  by  SESAM  is 

a  =  (\-  0)as  +  da  ads 

Where  8  is  the  effective  fractional  coverage  of  atomic  oxygen,  aads  is  the  accommodation  coefficient  of 
the  adsorbate  (atomic  oxygen),  typically  assumed  to  be  unity  for  complete  accommodation,  and  as  is 
calculated  with  Ks  of  2.4.  6  is  calculated  in  a  similar  manner  to  a  in  the  original  model  ECAM,  except  for 
a  few  minor  changes  which  are  not  discussed  here. 

The  modified  semi-empirical  satellite  accommodation  model  (SESAM  (MOD))  is  essentially  the  same 
model  as  SESAM  with  3.6  as  the  value  for  Ks  in  Goodman’s29  formula. 

5.4  Results 

The  developed  drag  coefficient  model  for  GRACE  was  used  to  compute  CD  for  the  duration  of  the 
selected  days.  Comparison  of  CD  computed  using  the  different  energy  accommodation  coefficient  models 
and  constant  energy  accommodation  coefficient  of  0.93  (used  by  Sutton)  is  performed.  Precision  attitude 
and  orbit  data  was  provided  by  Sutton.  CD  for  attitudes  that  are  within  the  bounds  used  in  model 
development  is  computed  using  linear  interpolation.  CD  for  attitudes  that  are  beyond  the  bounds  used  in 
model  development  is  computed  using  linear  extrapolation.  Precise  orbit  data  was  used  to  obtain 
atmospheric  properties  (Tinf  and  alpha)  from  NRLMSISE00. 

Figure  14  shows  the  CD  computed  for  GRACE  using  the  developed  model  for  July  19,  2005  and 
simulated  solar  max  conditions.  Simulation  of  solar  max  conditions  uses  a  10.7  cm  so-lar  radio  flux  index 
(FI 0.7)  value  of  260  and  a  geomagnetic  disturbance  index  (Ap)  value  of  21 .  The  percent  bias  in  CD 
computed  with  a  constant  energy  accommodation  coefficient  value  of  0.93  lie  in  the  range  of  -7%  to  5% 
with  the  mean  close  to  0%.  The  percent  bias  in  CD  computed  using  EACM  (MOD)  lie  in  the  range  of  - 
8%  to  3%  with  the  mean  close  to  0%.  The  percent  bias  in  CD  computed  SESAM  lie  in  the  range  of  -4% 
to  7%  with  the  mean  close  to  2%.  The  percent  bias  in  CD  computed  SESAM  (MOD)  lie  in  the  range  of  - 
7%  to  5%  with  the  mean  close  to  0%.  Since,  the  constant  value  of  0.93  for  energy  accommodation 
coefficient  simulates  solar  max  conditions,  the  mean  different  between  CD  computed  by  Sutton  and 
DSMC  using  different  energy  accommodation  coefficient  models  is  around  0%.  The  range  and  mean  of 
the  percent  bias  between  CD  and  density  computed  using  different  models  and  Sutton  will  henceforth  be 
denoted  as  (range,  mean).  The  ECAM  (MOD)  and  SESAM  (MOD)  show  identical  results  because  they 
use  the  same  Goodman’s  formula  with  different  formulations.  The  mean  bias  in  CD  computed  using  the 
SESAM  model  is  higher  than  the  other  models  because  it  uses  2.4  for  Ks  in  Goodman’s  formula. 
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Figure  14:  (top)  CD  computed  using  different  energy  accommodation  coefficient  models  for  July  19,  2005 
and  simulated  solar  max  conditions,  (bottom)  Percent  different  in  CD  computed  using  different  energy 
accommodation  coefficients  model  and  CD  computed  by  Sutton’s  using  a  flat  plate  model. 

Figure  15Figure  shows  the  CD  computed  for  GRACE  using  the  developed  model  for  July  19,  2005  and 
simulated  solar  min  conditions.  Simulation  of  solar  min  conditions  uses  a  F10.7  value  of  50  and  an  Ap 
value  of  2.  The  percent  bias  in  CD  computed  with  a  constant  energy  accommodation  coefficient  value  of 
0.93  lie  in  the  range  of  -18%  to  -5%  with  the  mean  close  to  -10%.  The  percent  bias  in  CD  computed  using 
EACM  (MOD)  lie  in  the  range  of  3%  to  17%  with  the  mean  close  to  12%.  The  percent  bias  in  CD 
computed  SESAM  lie  in  the  range  of  7%  to  19%  with  the  mean  close  to  15%.  The  percent  bias  in  CD 
computed  SESAM  (MOD)  lie  in  the  range  of  3%  to  17%  with  the  mean  close  to  12%.  The  ECAM 
(MOD)  and  SESAM  (MOD)  again  show  identical  results  with  the  mean  bias  in  CD  computed  using  the 
SESAM  and  Sutton  again  being  higher. 
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Figure  15:  (top)  CD  computed  using  different  energy  accommodation  coefficient  models  for  July  19,  2005 
and  simulated  solar  min  conditions,  (bottom)  Percent  different  in  CD  computed  using  different  energy 
accommodation  coefficients  model  and  CD  computed  by  Sutton’s  using  a  flat  plate  model. 


Figure  16  shows  the  density  computed  for  GRACE  by  scaling  accelerometer-derived  densities  from 
Sutton  using  appropriate  CD  and  area  ratios  for  July  19,  2005  and  simulated  solar  min  conditions.  The 
percent  bias  in  density  computed  with  a  constant  energy  accommodation  coefficient  value  of  0.93  lie  in 
the  range  of  2%  to  15%  with  the  mean  close  to  5%.  The  percent  bias  in  density  computed  using  EACM 
(MOD)  lie  in  the  range  of  -18%  to  -8%  with  the  mean  close  to  -15%.  The  percent  bias  in  density 
computed  SESAM  lie  in  the  range  of  -1 1%  to  2%  with  the  mean  close  to  -6%.  The  percent  bias  in  density 
computed  SESAM  (MOD)  lie  in  the  range  of  -18%  to  -11%  with  the  mean  close  to  -16%. 
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Figure  16:  (top)  Sutton’s  accelerometer-derived  densities  scaled  with  appropriate  CD  and  area  ratios  using 
different  energy  accommodation  coefficient  models  for  July  19,  2005  and  simulated  solar  min  conditions, 
(bottom)  Percent  different  in  density  calculated  using  different  energy  accommodation  coefficients 
models  and  accelerometer-derived  densities  by  Sutton. 


5  Conclusions 

Precision  orbit  ephemeris  (POE)  data  are  used  as  measurements  in  an  orbit  determination  process  to 
estimate  neutral  atmospheric  density  along  the  CHAMP  and  GRACE  orbits.  The  combination  of  a  CIRA 
1972  as  a  baseline  density  model  with  ballistic  coefficient  half-life  of  1.8  minutes  and  density  half-life  of 
1 80  minutes  was  found  to  be  the  overall  combination  for  estimating  POE  derived  densities  that  best 
matched  accelerometer  derived  densities  as  measured  by  cross  correlation  and  RMS. 

POE  derived  densities  were  found  for  the  entire  life  of  CHAMP  and  from  launch  through  2011  for 
GRACE.  The  density  data  exists  as  a  continuous  data  set  divided  into  one  week  files  for  ease  of  access. 
The  POE  derived  densities  were  shown  to  better  match,  in  terms  of  cross  correlation  and  RMS,  the 
accelerometer  derived  densities  than  the  HASDM,  NRLMSISE-00,  and  Jacchia  71  density  models.  The 
cross  correlation  between  other  density  sources  and  the  accelerometer  derived  densities  was  lower  when 
the  orbit  plane  of  the  satellite  was  near  the  terminator  because  the  high  frequency  variations  observed 
only  by  the  accelerometer  become  similar  in  magnitude  to  the  daylight/darkness  variations  during  an 
orbit. 

Several  times  during  the  lifetime  of  CHAMP  and  GRACE  the  satellites’  orbits  were  nearly  coplanar.  One 
of  those  time  periods  occurred  during  April  2005.  For  a  few  hours  during  April  3  and  again  on  April  5, 
the  GRACE  satellites  were  flying  nearly  directly  above  CHAMP.  A  comparison  was  made  of  the  density 
estimated  from  CHAMP  and  GRACE-A  during  these  two  days  using  accelerometer,  POE  density, 
HASDM,  and  the  Jacchia  1971  and  NRLMSISE-00  models.  NRLMSISE-00  and  Jacchia  1971 


overestimated  the  densities  for  both  CHAMP  and  GRACE  during  these  time  periods.  The  POE  and 
HASDM  densities  matched  the  overall  changes  of  the  accelerometer  quite  well  for  CHAMP,  but 
overestimated  the  sunlit  density  peaks  for  GRACE.  Also,  there  are  high  frequency  variations  that  are  only 
observed  by  the  accelerometers.  An  interesting  feature  observed  in  the  accelerometer  densities  is 
nocturnal  density  peaks  that  occur  during  both  days  that  correspond  to  the  Midnight  Density  Maximum 
phenomena.  These  peaks  are  either  smoothed  through  or  show  up  as  slight  humps  or  increases  in  density 
trough  magnitudes  for  the  POE  and  HASDM  densities,  but  are  modeled  by  NRLMSISE-00. 

Densities  were  estimated  for  CHAMP,  GRACE,  and  TerraSAR-X  for  September  26-27,  2007  and 
September  29-30,  2007.  The  POE  densities  were  able  to  accurately  observe  the  major  changes  seen  in  the 
CHAMP  and  GRACE  accelerometers,  but  not  the  high  frequency  variations.  Jacchia  1971  significantly 
overestimated  the  CHAMP  and  GRACE  densities  during  these  time  periods.  For  TerraSAR-X,  at  a  higher 
altitude,  the  POE  and  Jacchia  1971  densities  were  very  close.  This  is  either  because  Jacchia  1971  better 
models  the  density  at  higher  altitude  than  at  lower  altitude  or  because  the  density  variations  are  less 
observable  using  the  POE  density  technique  because  of  the  significantly  lower  density  magnitudes.  The 
strong  MDM  seen  in  the  CHAMP  and  GRACE  accelerometer  densities  on  April  3  and  April  5  of  2005 
were  observable  in  the  September  26-27  and  September  29-30  of  2007  CHAMP  accelerometer  densities 
as  well.  However,  they  are  not  observed  in  the  GRACE  accelerometer  data  or  in  any  of  the  POE  densities 
for  CHAMP,  GRACE,  or  TerraSAR-X. 

The  effects  of  limited  data  and  lower  quality  data  on  density  estimation  are  examined  to  see  the  suitability 
of  using  a  larger  set  of  satellites  to  estimate  density.  Discontinuous  data  sets  are  created  by  decimating 
the  original  POE  data  for  CHAMP  and  GRACE.  Sets  were  created  with  15  and  8  minutes  of  data  per 
orbit.  Noise  levels  from  10  cm  to  100  m  were  added  to  the  data  to  simulate  data  of  lower  accuracy.  The 
results  showed  decreases  in  cross  correlation  and  RMS  of  the  resulting  density  estimations  compared  to 
accelerometer  derived  densities.  However,  the  data  even  for  relatively  high  noise  levels  should  be  usable 
for  estimating  density.  This  means  many  more  potential  satellites  can  be  used  to  estimate  density 
providing  much  better  spatial  coverage  of  the  thermosphere. 

The  effects  of  high  frequency  variations  in  atmospheric  density  on  orbit  propagation  have  been  examined. 
The  high  frequency  variations  are  observed  in  densities  derived  from  the  accelerometers  onboard  the 
CHAMP  and  GRACE  satellites,  but  are  not  observed  in  densities  derived  from  precision  orbit  ephemeris 
(POE)  data,  or  the  Jacchia  71  or  HASDM  density  models.  Specific  high  frequency  variations  examined 
were  from  Traveling  Atmospheric  Disturbances,  polar  cusp  density  enhancements,  and  the  times  when 
the  orbit  plane  was  near  terminator  and  high  frequency  variations  are  of  the  same  order  as 
daylight/darkness  variations.  The  POE  derived  densities  were  superior  for  orbit  propagation  compared  to 
Jacchia  71  in  all  cases  and  superior  or  comparable  to  HASDM  in  all  cases.  The  errors  between  the  orbits 
propagated  with  POE  derived  densities  compared  to  orbits  propagated  with  accelerometer  derived 
densities  were  around  1  m  for  GRACE  and  ranged  from  4-22  m  for  CHAMP.  These  errors  are  not 
significant  for  most  orbit  applications.  Therefore,  the  high  frequency  variations  observed  only  by  the 
accelerometer  are  not  significant  for  most  orbit  analysis  applications. 

A  model  to  accurately  predict  the  drag  coefficient  along  the  orbit  of  the  GRACE  satellite  was  successfully 
developed  using  the  Direct  Simulation  Monte  Carlo  technique.  The  effect  and  importance  of  accurate 
geometry  modeling  on  drag  coefficient  has  been  examined  using  the  developed  model.  Drag  coefficients 
computed  with  the  model  using  different  energy  accommodation  coefficients  models  are  compared  with 
those  derived  by  Sutton  using  a  flat  plate  model  and  a  constant  energy  accommodation  coefficient  of  0.93. 
Drag  coefficients  computed  with  the  model  using  a  constant  accommodation  coefficient  compared  with 
Sutton  show  a  mean  bias  in  the  range  of  3-5%  with  time  varying  biases  going  as  high  as  15%  resulting 
from  inaccurate  modeling  of  the  geometry.  Drag  coefficients  computed  with  the  model  using  different 
energy  accommodation  co-efficient  models  show  a  mean  bias  in  the  range  of  5-10%  with  time  varying 


biases  going  as  high  as  19%.  Use  of  a  constant  energy  accommodation  coefficient  and  a  flat  plate  model 
results  in  drag  coefficients  that  are  consistently  and  significantly  under  predicted  during  solar  minimum 
conditions.  During  solar  maximum  conditions,  the  mean  bias  in  drag  coefficients  computed  with  the 
model  and  by  Sutton  is  very  close  to  zero;  however,  time  varying  biases  still  reach  as  high  as  6-7%. 

The  effect  of  biases  in  drag  coefficient  on  density  estimation  using  accelerometer  data  has  also  been 
examined.  Accurate  density  estimates  have  been  calculated  by  scaling  the  accelerometer-derived  densities 
from  Sutton  using  appropriate  drag  coefficient  and  area  ratios  and  different  energy  accommodation 
coefficients  models.  Densities  calculated  with  the  model  using  a  constant  accommodation  coefficient 
compared  with  Sutton  show  a  mean  bias  in  the  range  of  2-3%  with  time  varying  biases  going  as  high  as 
12%  resulting  from  inaccurate  modeling  of  the  geometry.  Densities  calculated  with  the  model  using 
different  energy  accommodation  coefficient  models  show  a  mean  bias  in  the  range  of  8-12%  with  time 
varying  biases  going  as  high  as  18%.  Use  of  a  constant  energy  accommodation  coefficient  and  a  flat  plate 
model  results  in  drag  coefficients  that  are  consistently  and  significantly  over  predicted  during  solar 
minimum  conditions.  During  solar  maximum  conditions,  the  mean  bias  in  densities  calculated  with  the 
model  and  by  Sutton  is  in  the  range  of  4-6%  with  time  varying  biases  still  reaching  as  high  as  11%. 
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